Analytic study of the three-urn model for separation of sand 
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We present an analytic study of tlie three-urn model for separation of sand. We solve analytically 
the master equation and the first-passage problem. We find that the stationary probability distribu- 
tion obeys the detailed balance and is governed by the free energy. We find that the characteristic 
lifetime of a cluster diverges algebraically with exponent 1/3 at the limit of stability. 
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I. INTRODUCTION 

A granular system consisting of macroscopic particles 
exhibits extremely rich phenomena, which has been stud- 
ied both experimentally and theoretically One of 
such interesting phenomena is the spatial separation of 
shaken sand. In the experiment by Schlichting and Nord- 
meier , granular particles are prepared in a box which is 
mounted on a shaker and separated into two equal parts 
by a wall. There is a slit on a wall through which parti- 
cles can move from one compartment to the other. Under 
a certain shaking condition, the granular particles simul- 
taneously separate into dense and dilute regions, which 
will not occur for gaseous particles, due to the dissipative 
nature of the macroscopic particle collision. 

The emergence of symmetry breaking in shaken sand 
was first explained by Eggers using a hydrodynamic ap- 
proach 3], and later by Lipowski and Droz using an urn 
model which is supposed to capture the essence of the 
experimental system In the urn model, N granular 
particles are distributed into L urns. And each particle 
can jump from one urn to another with the probability 
controlled by a parameter, called an effective tempera- 
ture, which depends explicitly on the density of particles 
in each urn. The dissipative nature of the particle col- 
lision is incorporated into the model with the density- 
dependent effective temperature. 

The urn model is simple enough to allow extensive 
numerical calculations 4, 5] as well as analytical stud- 
ies 0, ■ The two- urn {L = 2) model shows a rich struc- 
ture with symmetric, mixed, and asymmetric phases sep- 
arated with continuous and discontinuous transitions as 
well as the tricritical point. In the symmetric phase par- 
ticles are distributed equally in each urn, while the sym- 
metry is broken in the asymmetric phase. In the mixed 
phase, both the symmetric and asymmetric states are sta- 
ble. It was also found that the characteristic time which 
it takes to reach the symmetric state from an asymmetric 
state is given by the same free energy function which gov- 
erns the stationary probability distribution Coppex 
et al. also numerically investigated the three-urn model 
to find the absence of continuous transition They 
also obtained that the characteristic time at the limit of 
stability diverges algebraically as the number of particles 
increases (see also Ref. Q)- 



In this paper, we present the results of an analytic 
study to the master equation and the characteristic time 
scale in the three-urn model. We solve the master equa- 
tion in the thermodynamic limit to find that the solution 
shows a nature of deformed wave. We also obtain the sta- 
tionary probability distribution in a finite system, which 
interestingly obeys the detailed balance. Combining the 
knowledge of the characteristic scales of those distribu- 
tion and the mean-field flux equation, we obtain the scal- 
ing property of the characteristic time scale at the limit 
of stability. 

The paper is organized as follows. In Sec. ^] we briefly 
review the model and its master equation. In Sec. IIIII 
we present the analytic solution of the master equation 
in the thermodynamic limit and the analytic expression 
of the stationary probability distribution. In Sec. IIVI we 
investigate the scaling law of the characteristic time scale. 
Section [3 is devoted to conclusions and discussions. 



II. MODEL AND ITS MASTER EQUATION 

The model introduced by Coppex et al. is deflned as 
follows. N particles are distributed between three urns, 
and the number of particles in each urn is denoted as Ni , 
and N3 = N — Ni — N2, respectively. At each time of 
updates one of the N particles is randomly chosen. Let 
n be a fraction of the total number of particles in the urn 
which the selected particle belongs to. With probability 
exp(— j^^) the selected particle moves to a randomly 
chosen neighboring urn. T{n) — Tq + A(l — n) is the 
effective temperature of an urn with particles nN that 
measures the strength of fluctuations in the urn. For 
more detailed description of the model, we refer readers 
to Ref. 3. 

It is easy to derive the master equation for the proba- 
bility distribution p{Ni, N2,N3, t) that there are Ni par- 
ticles in urn i at time t 

p{Ni,N2,N3,t+l) = 

p{^^)l(piNi + i,N2-i,m,t) 

+p{Ni + l,N2,N3-l,t)) 
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p(7Vi,7V2,iV3,t), 



(1) 



where F{n) = nexp(— measures the flux of parti- 
cles leaving the given urn. Here we introduced for conve- 
nience the notation p(A^i, 7V2, ^3, = if iVi-|-iV2-l-iV3 ^ 
N or any A^^'s is either less than or greater than N . 

Let's denote the occupancies of the urns by rii — Ni/N . 
The time evolution of the averaged particle occupancies 
is governed by the equations 



where 
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M)t+i = {nr)t + -^{^i{n))t 



(2) 



/* - Y.N„N,,N,{---)PiNi^N2,N^.t), and 

i X]fc=i ^(^fc) ~ ^F{ni). The unit of time may 
be chosen so that there is a single update per a particle on 
average. Scaling the time by iV, expanding Eq. Q with 
respect to , and using the mean- field approximation in 
evaluating the average, we get 



(3) 



where n.^^t) = {ni)t. 

Detailed analysis on the existence of the stable station- 
ary solutions of Eq. ^ was done by Coppex et al. 0. 
We here display their phase diagram in Fig.^to make our 
paper as self-contained as possible. The stable symmet- 
ric solution (ni = n2 =1^3) exists in region I, III, and IV 
while the stable asymmetric solution (ni > n2 ~ na) ex- 
ists in region II, III and IV. Note that Eq. Q is obtained 
with mean- field approximation that {J-{n)) = !F{{n)). 
However, as will be shown later, the relation becomes ex- 
act for the delta-peaked initial probability distribution, 
so is the phase diagram. 



III. THE SOLUTION OF THE MASTER 
EQUATION 

We are mainly interested in investigating the proper- 
ties of the infinite system. Consider the thermodynamic 
limit N ^ 00 with the occupancies of the urns being 
fixed. Representing the probability distribution by in- 
stead of A^i, scaling the time by TV, expanding Eq. 
with respect to and keeping the terms up to the first 
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FIG. 1: Phase diagram of the three-urn model The sym- 
metric solution vanishes continuously on the solid line while 
the asymmetric one disappears discontinuously on the dotted 
line. The transition of the behavior of the stationary prob- 
ability distribution is denoted by the dashed line. For the 
inset, see Sec. II VI 



order, we arrive at the partial differential equation 



d 



dt 



p{n, t) 



d d 
— [Ti{n)p{n, t)] + — [T2{n)p{ri, t)] - . 

(4) 

Here we would like to remind that 77.3 = 1 — ni — n2 so 
that the independent variables are ni,n2. 

Note that Eq. ^ may be interpreted as the continuity 
equation for the probability with velocity {Ti^ T-2)- Since 
the velocity is already given as a function of n, the solu- 
tion to Eq. Q can be formally found as follows. Time 
evolution of a point located at r = (a;, y, 1 — a; — y) at 
t = is determined by the Eq. We denote it by 

R{t;r). We show its typical trajectories in Fig. El In 
the symmetric phase (region I in Fig. ^ there is only 
the stable fixed point corresponding to the symmetric 
state so that every trajectory flows to that point. This is 
shown in Fig. |2Ia). Figure|2Ib) shows a typical behavior 
for the asymmetric phase (region II in Fig. In this 
case, there are three stable fixed points corresponding to 
asymmetric states and one unstable fixed points corre- 
sponding to the symmetric state as well as three saddle 
points. The trajectories are separated by the separatri- 
ces denoted by dashed lines. Finally Fig. [2Ic) shows a 
typical behavior for the mixed phase (region III and IV 
in Fig. ^ . There are one stable symmetric fixed point 
and three asymmetric fixed points as well as three saddle 
points. The trajectories flowing toward the stable fixed 
points are separated by the separatrices. 

Since the continuity equation implies that the proba- 
bility at point r evolves according to i?(i; r), we obtain a 
formal solution of Eq. Q 



R{t-T) 



(5) 



3 




FIG. 2: Typical trajectories of Eq. (|^ . We plotted the trajec- 
tories in the diagonal plane of a unit cube so that each variable 
rii are treated in the same way. (a) A = 0.5, To = 0.5 (sym- 
metric phase) (b) A = 0.25, To = 0.05 (asymmetric phase) 
(c) A = 1.0, To — 0.1 (mixed phase). The stable, the unsta- 
ble, and the saddle point are represented with a filled circle, 
a filled square, and a cross, respectively. 



It implies that the initial probability distribution in the 
basin of attraction associated with a stable fixed point 
will eventually accumulate at that point. As a conse- 
quence, in the long time limit f — > oo the probability 
distribution becomes a sum of delta peaks at the stable 



fixed points denoted by nf, 

p{n, oo) = ^ piS^"^^ {n- Hi) , 



(6) 



where pi are the sum of the initial probabilities in the 
basin of attraction associated with n^. 

Now we will consider another limit in the master equa- 
tion (^), namely take the long time limit t oo before 
we take the limit TV — > cx). That is, we are looking for 
the stationary probability distribution for a finite system. 
Though this limit may not properly reflect the properties 
of the infinite system, we expect it will reveal interest- 
ing properties of the system concerning the characteristic 
times (See @ in two-urn case.). 

Let's first take the long time limit of t — > oo in Eq. 
In this limit we may drop off the time dependence in the 
probability distribution, which now takes the form 



^ JVi + 1 p(iVi-H,jV2-l,jV3) _ 
^ N ' p{Ni,N2,N3) 
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N ' p{N,,N2,N3) 



0(7) 



In contrast to the corresponding equation (13) in the 
two- urn model , it does not show a simple tridiagonal 
structure. However, we can show that the stationary 
probability distribution determined by Eq. Q obeys the 
detailed balance 



■ N 



,N. 



)p{N[,N^,N'^) = F{^)p{N,,N2,N:,) , (8) 



where {i,3,k} = {1,2,3} and N[ ^ Ni + 1, A^j = - 
l,iV^ = Nk- We also would like to note that this kind 
of relation is obeyed in two-urn model. Equation Q 
essentially tells us that probability distribution with Nj 
is given by that with Nj — 1, which is in turn given by 
that with Nj — 2, and so on. Repeatedly using Eq. {Tj) 
with i as 1 or 2 and i = 3, we obtain 



piN,,N2,N3) 



N 



-p{0,0,N), 



(9) 

where p(0, 0, N) appears as an overall factor to normalize 
the probabilities so that we get 



p(0,0,iV) = 



1+ 



N 



N~Ni 



N 
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Now let's take the limit N 



^ = 712, and scaling the probability distribution by N'^, 
the stationary probability distribution for large N now 
becomes 



oo. With ^ = m, 



pin) 



(10) 



where 



/ dt[lnF(l-i)] 

dt [\nF{t)] - dt [\nF{t)] . (11) 



Gin) 



We will call G(n) the (negative) free energy function. 

In the limit TV — > oo, the main contribution to the 
stationary probability distribution comes only from the 
maximum of G(-), and it becomes the sum of delta peaks. 
The maximum of Gin) occurs when 



d 
dni 



Gin) 



d 
dn2 



Gin) = 



Fim) ^ Fin2) ^ Fins) . 



(12) 



(13) 



This condition is equivalent to the stationary condition 
of the flux equation JHJ. Note that in region II only 
the three stable asymmetric solutions having the same 
maximum exist, while in region I only the symmetric 
solution is stable. Therefore the stationary probabil- 
ity distribution has the triple peaks in region II and 
only the central peak in region I. In region III and IV 
both the symmetric and the asymmetric solutions are 
stable so that the maximum of G(n) should be deter- 
mined by comparing its values at the stable fixed points. 
The crossover of the maximum point occurs when both 
values coincide. This implies that the transition be- 
tween the triple peaks and the central single peak in the 
probability distribution is determined by the condition 
AG = G(na) - G(l/3, 1/3, 1/3) = 0, where Ha is one of 
the stable asymmetric fixed points. This condition yields 
a line that separates two regions III and IV. 



IV. CHARACTERISTIC TIME SCALE 

One of the interesting phenomena in sands separated 
by compartments is a sudden collapse of a granular clus- 
ter [^, ^] . It is observed experimentally that a granular 
cluster with majority of sand grains in a single compart- 
ment remains stable for a long time until it collapses 
abruptly and diffuses over all compartments j^. The 
urn model proposed in Ref. Q displays the same phe- 
nomenon: In region I a granular cluster is stable up to 
a time scale t, after which sand particles are distributed 
uniformly over all boxes. Approaching the phase bound- 
ary I-IV, the characteristic time scale diverges. At the 



phase boundary, it is found numerically that r scales as 
T with z ~ 0.32 jSj]. In this section, the scaling 

law for the characteristic time r at and near the phase 
boundary is derived analytically. 

The master equation in Eq. Q in the large N limit 
does not contain a diffusive term. It implies that 
a delta-peaked probability distribution remains delta- 
peaked during time evolution, which enabled us to write 
down the formal solution in Eq. |SJ|. The dispersion-less 
nature also guarantees that the mean field approximation 
in Eq. ijjjl for the occupancy becomes exact as long as 
the initial values of rii's are prescribed. Hence we can 
study the dynamics of the granular cluster using Eq. Q 
with the initial condition rii = n2 = and n^ — 1. 

With ni = n2 = n and n^ = 1 — 2n, one can rewrite 
Eq. Q as 



V{n) 



2n) - Fin)} 



(14) 



with the initial condition n(i = 0) = 0. The cluster 
dynamics is then determined from the property of the 
flow function Vin). In the inset of Fig. ^ we show the 
plots of the flow function at different values of Tq with 
A — 0.3 flxed. At A = 0.3, the critical point separating 
the region I and VI is given by Tq = Tqc = 0.1698- • •. 
For To > Tqc, n = 1/3 is the unique stationary state 
solution; n grows from zero to 1/3 to reach the symmetric 
state with rii = n2 = n^ — 1/3. Note that there exists 
a local minimum in Vin) at < no < 1/3 where the 
flow velocity is minimum. Hence the value of n remains 
at the intermediate value n n^ for a long time, and 
then converges to n = 1/3 quickly. It turns out that the 
sudden collapse of a granular cluster H, 11] is due to the 
local minimum in Vin). The characteristic time scale or 
the life time r is given by 



dn 



(15) 



As Tq Tq^, the minimum flow velocity Vin^) decreases 
and hence r increases. The asymmetric solution with 
ni — n2 < 7i3 emerges at T — Tqc when Vino) — and 
T diverges. 

The scaling law for the characteristic time scale is 
determined from the analytic property of Vin) near 
Tq ~ Tqc and n ~ tiq. One can expand as Vin) ~ 
a{n ~ nQ)'^ + b{TQ — Tqc) , where a and b are constants of or- 
der one and all higher order terms in (Tq ~Tqc) are irrele- 
vant. Inserting it into Eq. H15|l . one obtains that the char- 
acteristic time scale diverges as r ~ (Tq — Tqc)^^^^ (see 
also Ref. Q). At To = Tqc, the granular system ap- 
proaches the asymmetric state algebraically in time as 
\n — no\ ~ whose life time r diverges. That is to say, 
the asymmetric state is stable in the N —>■ oo limit. 

However, for finite N, the system could evolve into the 
symmetric state at To — Tqc with the help of statisti- 
cal fluctuations in n. Indeed, the probability distribu- 
tion p has a flnite dispersion for flnite N, even though 



it becomes a sum of delta peaks for infinite N. In or- 
der to estimate the order of magnitude of fluctuations 
in n, we investigate the stationary probabihty distribu- 
tion in Eq. (|10|l at Tq = Tqc in detail. Near the asym- 
metric state with ni — n2 = n and 713 = 1 — 2n, the 
(negative) free energy function in Eq. Ullfl is written as 
G(n) = /p" dtlnF{l - t) - 2 £ dt\nF{t). Now we ex- 
pand it near n — uq to yield G{n) — G(rto) + G'{no){n — 
710 ) + ^G"{na){n ~ no f + ^G"'{no)in - no)^ + • • •• After 
straightforward calculations, one can show that 

G'K) - 21n(i^(l-2no)/F(no)) (16) 

G"M - ■J—,v'M (17) 

G"'K) = t;^T^"K) . (18) 
F[no) 

Note that G'(no) = G"{no) = since V{na) = V'{no) = 
0. Therefore, the (negative) free energy function can be 
approximated as G{n) ~ 0(7^— no)^ with a constant c, and 
hence the probability distribution as Ps{n) ~ gcN{n-no) 
up to a normalization constant. 

The stationary probability distribution suggests that 
the magnitude of the fluctuation in n near uq is of order 
5n N~^/'^. With the help of the fluctuation, the gran- 
ular system could get away from the asymmetric state 
when \n{t) — 7^o| ~ Sn, and flow into the symmetric 
state. Therefore the characteristic time scale is given 

by r ~ /o"""^" dt/V{t), which resuhs in r ~ with the 
dynamic exponent 

z = l/3. (19) 
Our analytic result confirms the numerical result z 2± 0.32 
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reported in Ref. j^. 



V. CONCLUSIONS 

We analytically investigate the three-urn model intro- 
duced by Coppex et al. [3 We formally solve the master 
equation of the model in the thermodynamic limit and 
find how the probability distribution evolves. In the long 
time limit, the probability distribution becomes delta 
peaks only at the stable fixed points. The strength of 
a delta peak is equal to the sum of initial probabilities in 
the basin of attraction associated with that fixed point. 

We solve exactly the stationary probability distribu- 
tion where we take the long time limit before we take 
thermodynamic limit. Wc find the distribution obeys 
the detailed balance. Regardless of the initial probability 
distribution it shows triple peaks or a single central peak 
depending on the parameters of the system. The final 
formula of the stationary probability distribution resem- 
bles that of the equilibrium systems, where the transition 
from the triple peaks to the single peak is determined by 
the condition that free energies of two phases become 
equal. 

We also obtain the exact scaling law for the character- 
istic time scale r which it takes to reach the symmetric 
state from an asymmetric state near and at the phase 
boundary I- VI. In the symmetric phase (region I), the 
granular cluster is unstable and t is finite. It grows as one 
approach the phase boundary as t ~ (Tq — Tqc)^^^^. At 
To — Tqc, the characteristic time diverges algebraically 
as T ~ with the dynamic exponent z = 1/3. 
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